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Abstract 

The statistical properties of an ecosystem composed of species in- 
teracting via pairwise, random interactions and deterministic, concen- 
tration limiting self-interactions are studied analytically with tools of 
equilibrium statistical mechanics of disordered systems. Emphasis is 
given to the effects of externally induced extinction of a fixed fraction 
of species at the outset of the evolutionary process. The manner the 
ecosystem copes with the initial extinction event depends on the de- 
gree of competition among the species as well as on the strength of 
that event. For instance, in the regime of high competition the ecosys- 
tem diversity, given by the fraction of surviving species, is practically 
insensitive to the strength of the initial extinction provided it is not 
too large, while in the less competitive regime the diversity decreases 
linearly with the size of the event. In the case of large extinction events 
we find that no further biotic extinctions take place and, furthermore, 
that rare species become very unlikely to be found in the ecosystem at 
equilibrium. In addition, we show that the reciprocal of the Edwards- 
Anderson order parameter yields a good measure of the diversity of 
the model ecosystem. 



PACS: 87.10+e, 87.90.+y, 89.90+n 



1 Introduction 



Extinction seems to be the final outcome of the evolution of species. In fact, as 
species survive about 10 million years in the average, nearly all species that have 
ever existed are extinct and only a very small fraction of them have left their impres- 
sions in the fossil record [El . The causes of the mass extinction events is currently 
a matter of dispute as there are two main types of explanations || . The more tra- 
ditional one asserts that extinction is caused by external stresses as, for instance, 
major climate changes and asteroids impacts. This point of view is supported by 
some evidences such as the unusual quantity of iridium and other noble metals in 
the rocks that marked the boundary between the Cretaceous and Tertiary periods, 
when the era of the dinosaurs was replaced by the era of the mammals [Q . Since 
iridium is more common in asteroids than in the Earth's crust, this finding can be 
viewed as evidence for an asteroid impact. The alternative explanation asserts that 
extinctions are caused by the interactions between the species in the ecosystems. 
In particular, Paine || has shown that species richness sometimes can be increased 
by the predator-mediated coexistence, and the removal of predators can lead to 
additional species extinctions. Some recent studies indicate that food webs with 
many species or high connectivity are more likely to lose species as a consequence 
of the extinction of a single species when compared with more simple food webs 
H §. Although this kind of argument seems well suited to explain the so-called 
background extinctions, it certainly needs some new ingredients to explain mass 
extinctions as well. In fact, the missing ingredient seems to be the self-organized 
criticality concept || , which in this context is best illustrated by the popular Bak- 
Sneppen model ||. According to this model, the fitness of each species is affected 
by the other species to which it is coupled in the ecosystem so that large events 
in the evolutionary history may be thought of as large coevolutionary avalanches 
caused by the intrinsic dynamics of the model. In this model, the distribution of 
the extinction sizes follows a power law, which is a valid candidate for fitting to the 
experimental data. 

Although we recognize that evolution and hence extinction are, as pictured by 
the models mentioned above, essentially dynamical phenomena, in this work we 
study these phenomena within the equilibrium statistical mechanics framework of 
the random replicator model for species coevolution jnj |ll[ |l2|, Deterministic 
replicator models are commonly used to describe the evolution of self-reproducing 
entities in a variety of fields such as game theory, prebiotic evolution and sociobiol- 



ogy |L4, 15 . The random replicator model introduced by Diederich and Opper [[L0| 
attempts to model the uncertainties and the overwhelming complexity of the in- 
terspecies interactions in biological ecosystems by assuming that those interactions 
are random. However, it is also assumed that the dynamics is such that a fitness 
functional (Lyapunov function) is maximized so that the only stationary states are 
fixed points. In fact, the existence of such a functional leads to a replicator equa- 
tion with symmetric interspecies interactions JT5[ which is a severe assumption from 
the biological standpoint. It allows, however, full use of tools of the equilibrium 
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statistical mechanics to study analytically the average properties of the equilibrium 
states of this kind of disordered ecosystems. 

An interesting result of the random replicator model is that in the equilibrium 
state a fraction of the species is extinct . The mechanism of extinction is clearly 
outcompctition and, in the absence of any cooperation pressure, only the pair of 
species with the largest reinforcing interactions will thrive. In this contribution we 
study the effects of random elimination of a fixed fraction of the species at the out- 
set of the evolutionary process, giving emphasis to the distribution of the remaining 
species concentrations in the equilibrium state. There are several interesting issues 
that can be addressed in this framework. For instance, what is the equilibrium situ- 
ation when the fraction of species eliminated at the beginning is already larger than 
the fraction that would be extinct naturally due to outcompctition? Furthermore, 
how does the ecosystem cope with large initial extinction events? In this paper we 
give clear-cut analytic answers to these questions which are partly corroborated by 
numerical simulations of the model ecosystem. 

The remainder of the paper is organized as follows. In Sec. we introduce the 
model and discuss the ecological interpretation of the control parameters. The equi- 
librium properties of the model are derived within the replica-symmetric framework 
in Sec. [| and the use of the reciprocal of the Edwards- Anderson order parameter 
as a measure for the diversity of the ecosystem is suggested. In Sec. |] we calculate 
the distribution of probability of the concentrations of a given species, allowing 
thus the explicit calculation of the ecosystem diversity as the fraction of surviving 
species at equilibrium. Finally, in Sec. [| we present some concluding remarks. 

2 Model 

We consider an infinite population (ecosystem) composed of individuals belonging 
to N different species whose fitness Ti (i = 1,...,N) are the derivatives T{ = 
dT jdxi of the fitness functional T defined as 

- T = H (x) = u bix\ + ^2 Jij biXibjXj (1) 

i i<j 

where Xi € [0, oo) is the fraction of species i and bi is a quenched random variable 
that takes the values and 1 with probabilities a and 1 — a, respectively. Hence Na 
randomly chosen species are eliminated at the outset in the average and so hence- 
forth we will refer to a 6 [0, 1] as the dilution parameter. An effective competition 
among the species is enforced by requiring that the concentrations of the surviving 
species satisfy the constraint 

N 

hxi = Q N, (2) 

i=l 

where Qo is an arbitrary positive constant which gives the scale of the concen- 
trations Xi. The coupling strengths Jy between species i and j are statistically 
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independent quenched random variables with a Gaussian distribution 



(JijfN 



(3) 



so that Jij < corresponds to pairs of cooperating species while > to pairs of 
competing species. The self-interaction parameter u > acts as a global coopera- 
tion pressure limiting the growth of any single species, and it is crucial to guarantee 
the existence of a nontrivial thermodynamic limit, N — > oo. In fact, for large u 
the minimum of H. corresponds to a homogeneous ecosystem where the surviving 
species have concentrations Xi/Qo = 1/(1 — a) Vi. The positive self-interactions 
means that individuals of a same species compete against themselves, which is 
quite reasonable as they certainly share the same resources (ecological niche). 

The time evolution of the species concentrations is given by the replicator equa- 
tion 



dxi 
~~dt 



&H(x) 

dxi 



1 ^ dH(x) 

k 



dx k 



Vi 



(4) 



which minimizes 7i(x) while keeping the term . biXi constant during the evolution. 
Hence the fixed points of this equation are the minima of 7i(x) and in the following 
we use the replica formalism to study analytically the statistical properties of these 
minima. 



3 The replica approach 

Following the standard prescription of performing quenched averages on extensive 
quantities only [l6| , we define the average free-energy density / as 

-^/ = Jim l(lnZ) (5) 

where 

z = I Il dx * s (q° n - 12 biX n \ Q N - 12 < l - e_/3W(x) ( 6 ) 

is the partition function and /3 = 1/T is the inverse temperature. Taking the limit 
T -> in Eq. (g) ensures that only the states that minimize TL (x) will contribute 
to Z. We impose the additional constraint 

-\h) = QN (7) 

i 

to avoid divergences when carrying out the integrals over Xi. Here (. . .) stands for 
the average over the coupling strengths Jy as well as over the auxiliary variables 
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bi. As usual, the evaluation of the quenched average in Eq. (|^) can be carried out 
through the replica method: using the identity 

(hxZ) = lim -hx{Z n ) (8) 

n—*Q n 

we first calculate (Z n ) for integer n, i.e., Z" = Y["l=iZ p , and then analytically 
continue to n = 0. The final result is 



Pf = 



i-o cxti i { E p p p p -yE^+T^ {pP)2 + E ^ + E w p 

P P P P<6 p 

E(« p5 ) 2 + E + E Pb ln G o(b,p p -, q ps , R p , Q p )\ (9) 

R n h—C\ s 



2 



p<8 p b=0 

where Pq = a, Pi = 1 — a, and 



• /o P !> p p-A 

-6 ^ A'a^ - (1 - 6) Q P x p y 



pS x p x s 



(10) 



We note that while we have calculated the average over the couplings Jy explicitly, 
we have used the self-averaging property jj J2i hiGo(&i) = J2b ^ b ^ nG o(b) to elimi- 
nate the site dependence of the hi variables. The relevant physical order parameters 
are 



i 
i 

which measure the overlap between a pair of different equilibrium states x p and 
x" 5 , and the overlap of an equilibrium state x p with itself, respectively. Here, (. . .)y 
stands for a thermal average taken with the Gibbs probability distribution 

W (x) = i 6 (^Q N - E ^) 5 (® N - E^ 1 - b ^ ex P H 87 * Wl • ( 13 ) 

To proceed further we assume that the saddle-point parameters are symmetric 
under permutations of the replica indices, i.e., p p — p, p p — p, q pS = q, q pS = 
q, R p — R and Q p = Q. With this prescription the evaluation of Eq. (||) is 
straightforward yielding the following replica-symmetric free energy density 

-Pf = -^-pQoR + a + a\n{9-) + l —^\n 



2 a 2 \2f3(2u~y) 



5 



£>2 , roo 

+0(1 - a) TIL + (l-a) Dzlneric 
2(2u-y) J-oo 



y/]3(R + Zy /g) 
y/2(2u-y) 



(14) 



where y = f3(p—q) and Dz = dz exp(—z 2 /2)/^/2~rr is the Gaussian measure. Already 
at this stage we can see that the concentration of species eliminated at the outset, 
given by the parameter Q, decouples from the other physical parameters and hence 
does not have any effect upon them. In the zero-temperature limit the saddle-point 
equations df/dq = 0, df/dy = and df/dR = are given by 

A = 2^(u-y), (15) 
2y(2u - y) = (1 - a) erfc (-A/V2) , (16) 

exp (-A 2 /2) - (2u - yf - (2u - y) (A 2 + l) y. (17) 



and 
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We note that the parameter associated to the concentration of surviving species 
Qo appears only as a scale of q and so henceforth we will set Qq = 1 without loss 
of generality. In the replica-symmetric framework the Edwards-Anderson order 
parameter q is defined by 

« = (^I>>tV ( 18 ) 

If the concentrations Xi were normalized to 1 rather than to N then q would give 
the probability that two randomly selected individuals are of the same species, a 
quantity known as Simpson's index ]l7j ]. Nevertheless, we can still give a simple 
physical interpretation to q. For instance, values of q of order of 1 indicate the 
coexistence of a macroscopic number of species (i.e., Xi s» 1 for an extensive number 
of species), while large values of q signalize the dominance of a few species only (i.e., 
Xi w N for a finite number of species). Of course, this interpretation is equivalent 
to that given above for Simpson's index, and so we can view 1/q as a measure 
of the diversity of the ecosystem. In Fig. [I] we present 1/g as a function of the 
dilution parameter a for several values of the cooperation pressure u. The results 
of the numerical solution of the replicator equation, Eq. (S) , for N = 500 are also 
presented. Each data point is the average over 100 realizations of the matrix of 
coupling strengths, starting with an uniform distribution of concentrations. Since 
the labeling of the species is arbitrary we can set 6^ = for i < aN and 6j = 1 
otherwise, without loss of generality. In addition, we choose a such that aN is 
integer for simplicity. In agreement with our interpretation, for a close to 1 wc 
can observe the vanishing of 1/q, which characterizes an ecosystem composed of 
a few species only. For small values of u the analytical results show the existence 
of a maximum of diversity for a nonzero value of the dilution parameter (see the 
inset in Fig. |l|); the numerical results however do not corroborate this finding. This 
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Figure 1: The diversity 1/q as a funcion of the dilution parameter a for 
(top to bottom) u = 1.3 (O).0.8 (v),0.6 (O),0.4 (□) and 0.3 (A). The 
symbols are the results of the numerical solution of the replicator equation. 
The inset highlights the region of the diversity maximum. 



discrepancy can be explained by the instability of the replica-symmetric solution. 
In fact, carrying out the standard local stability analysis we find that this 
solution is locally stable wherever the condition 

A = -l + — _L- -erfc ( -- ^= ) <0 (19) 
2(2u-y) 2 \ V2J 

is satisfied. Figure || shows the regions in the plane (a, u) where the replica- 
symmetric solution is stable. In particular, we find that for a — this solution 
is stable for u > l/v2 while for a = 1 it is stable for u > 1/2. Hence, the maxima 
observed in Fig. |l| are indeed artifacts of the replica-symmetry framework. Nev- 
ertheless, the agreement between the analytical and numerical results is already 
excelent for u > 0.6. The rather puzzling independence of the diversity on the 
dilution parameter for small u has a simple explanation as will be seen in the next 
section. 



4 Discussion 

Although the interpretation of the reciprocal of the Edwards-Anderson order pa- 
rameter as the ecosystem diversity yields some information on the distribution of 
species at equilibrium, a better understanding is achieved by calculating explicitly 
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Figure 2: Almeida-Thouless line separating the regions of stability (A < 0) 
and instability (A > 0) of the replica-symmetric solution. 

the probability distribution that the concentration of one of the (1 — a)N remaining 
species, say Xk, assumes the value x, defined by 



Vk (x) = lim 



Y[ d Xj b k 6 (x k - x) W (x) 



(20) 



with VV (x) given by Eq. (jl3|). As all non- vanishing species concentrations arc 
equivalent we can write Vk (x) = V (x) Vfe. Hence to evaluate Eq. (E(J) we introduce 
the auxiliary energy 



Haux (x) = H (x) 



h]>2 b k 5 (x k - x) 



(21) 



so that 



V(x) 



lim 



N/3 



(22) 



h=0 



/3-+DO Nfj dh 

where Z aux is the partition function (^) with 7i replaced by Ti. aux . Using Eq. 
( p2[ ) the calculations needed to evaluate V (x) become analogous to those used in 
the evaluation of the free-energy density (|l4|). In addition, to handle a possible 
singularity in the limit (3 — > oo it is more convenient to deal with the cumulative 
distribution function C (x) — J dx'V(x'). Carrying out the calculations within 
the replica-symmetric framework we obtain 



C(x) = {l-a)\l 



1 



-erfc 



1 / x(2u - 
V2 V 



y) 



A 



(23) 



where q, y and A are given by the saddle-point equations (p"5|)-(p"7|). In Fig. || we 
show C(x) for u = 0.8 and several values of a. The first point to note is that 
linij^oo C(x) = 1 — a yields the fraction of surviving species at the outset, as 
expected. In addition, a nonzero value of C(0) indicates that the probability dis- 
tribution V(x) has a delta peak at x = and so C(0) actually yields the fraction 
of the species that survived the initial externally induced extinction event but that 
were extinct later on due to outcompetition. In the regime of large dilution, say 
a > 0.8 in Fig. |^, the cumulative distribution is very small and practically constant 
for small concentrations, indicating that no further extinctions have taken place 
and, furthermore, that rare species are very unlikely to be found in the ecosystem 
at equilibrium. We note that the numerical simulations yield results practically in- 
distinguishable from the analytical ones. The rough independence of the diversity 
1/q on the dilution parameter a observed in Fig. [l] for small u is easily understood 
with the aid of the cumulated distributions. In fact, a direct measure of the ecosys- 
tem diversity is given by the fraction of surviving species 1 — a — C(0), which is 
shown in Fig. || as function of a. (We recall that a is the fraction of species that 
were extinct at the outset due to some external stress and C(0) is the fraction that 
died out due to outcompetition.) The remarkable similarity between these figures 
corroborates our interpretation of 1/q as a measure of the diversity. Clearly, the 
diversity is insensitive to variations of a whenever the fraction of extinct species in 
the undisturbed ecosystem (i.e. C(0) calculated at a = 0) is already considerably 
larger than a, so that the species eliminated at the outset would probably be extinct 
later on anyway. 



5 Conclusion 

Although the dynamics of the random replicator model may not look very appeal- 
ing, in the sense that it always leads to fixed points, the frustration caused by the 
competition between the concentration limiting self-interactions (u > 0) and the 
tendency to unlimited growth of pairs of strongly cooperative (Jy < 0) species 
results in a highly nontrivial equilibrium, characterized by many meta-stable states 
[|To| and a phase of replica symmetry breaking [pT[ . Of course, these very features 
make some aspects of the dynamics (e.g., slow relaxation and hysteresis effects) 
nontrivial as well. The wealth of ecologically relevant issues that can be addressed 
within this equilibrium framework can be appreciated, for instance, in the case of 
high-order interactions among the species where it has been reported the emergence 
of a threshold value which gives a lower bound to the concentration of the surviv- 
ing species, preventing then the existence of rare (low concentration) species in the 
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An important outcome of the equilibrium analysis of the random replicator 
model is the finding that in order to reduce the degree of frustration a fraction 
of the species dies out [ [I0[. This type of extinction has clearly a biotic cause, 
namely, outcompetition |19|. In this paper we study how the model ecosystem 
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Figure 3: Cumulative distribution of the concentration of the initially sur- 
viving species in equilibrium for u = 0.8 and (top to bottom) a = 0, 0.4, 
0.6, and 0.8. The dashed curves are the results of the numerical solution of 
the replicator equation. 

copes with abiotic or externally induced extinction, in which a fraction of randomly 
chosen species is eliminated at the beginning of the coevolutionary process. We find 
that in the regime of high competition (small it) the ecosystem diversity, i.e., the 
fraction of surviving species is practically insensitive to the strength a of the initial 
extinction provided it is not too large, while in the less competitive regime (large u) 
the diversity decreases linearly with increasing a. In the case of a large extinction 
event we find that no further (biotic) extinctions take place and, furthermore, that 
rare species become very unlikely to be found in the ecosystem at equilibrium. This 
is distinct from the result mentioned above for the case of high-order interactions 
where the probability of finding rare species in the ecosystem is strictly null |i~3| ] . 

An interesting by-product of our investigation is the finding that the recipro- 
cal of the Edwards- Anderson order parameter (i.e., the replica-symmetric overlap 
between two equilibrium states) serves as an easy-to-calculate measure of the di- 
versity of the model ecosystem. This opens the exciting possibility of interpreting 
the different hierarquical levels of the overlap order parameter in the full replica 
symmetry breaking scheme jl(| as different levels of a phylogenetic tree that gives 
the relations of dependence (viewed as ancestrality) among the species. 
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